# Install required libraries
# pip install numpy pandas matplotlib seaborn pyDOE scikit-learn catboost lightgbm xgboost

# Importing basic libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
import logging

# Suppress warnings for cleaner output
warnings.filterwarnings("ignore")

# Set logging level for LightGBM to ERROR to reduce verbosity
logging.getLogger('lightgbm').setLevel(logging.ERROR)

# Initialize a dictionary to store model performance metrics
modeldict = {
    'name': [],
    'MSE': [],
    'Rsquared': []
}

def addtodf(name, MSE, Rsquared):
    """
    Adds a new entry to the model performance dictionary.

    Parameters:
    - name (str): Name of the model.
    - MSE (float): Mean Squared Error of the model.
    - Rsquared (float): R-squared value of the model.
    """
    modeldict['name'].append(name)
    modeldict['MSE'].append(MSE)
    modeldict['Rsquared'].append(Rsquared)

# Using pandas to upload our dataset
df = pd.read_csv('data.csv', encoding='ISO-8859-1')

# Display the first few rows of the dataset
df

# Drop rows with missing values
df = df.dropna()
df

# Checking the structure and data types of the dataset
df.info()

<class 'pandas.core.frame.DataFrame'>
Index: 142 entries, 0 to 141
Data columns (total 32 columns):
 #   Column                                     Non-Null Count  Dtype  
---  ------                                     --------------  -----  
 0   Catalyst Used                              142 non-null    object 
 1   Interaction Type                           142 non-null    object 
 2   Support Interaction Level                  142 non-null    object 
 3   Support Type                               142 non-null    object 
 4   Catalyst Molecular Weight (g/mol)          142 non-null    object 
 5   Support Material Molecular Weight (g/mol)  142 non-null    float64
 6   Catalyst Electronegativity                 142 non-null    float64
 7   Support Material Electronegativity         142 non-null    float64
 8   Oxidative Level                            142 non-null    float64
 9   Reductive Level                            142 non-null    float64
 10  Surface Area (m_/g) min                    142 non-null    float64
 11  Surface Area (m_/g) max                    142 non-null    float64
 12  Pore Volume (cm_/g) min                    142 non-null    float64
 13  Pore Volume (cm_/g) max                    142 non-null    float64
 14  Optimal Temperature (K) min                142 non-null    float64
 15  Optimal Temperature (K) max                142 non-null    float64
 16  Particle Size (nm) min                     142 non-null    float64
 17  Particle Size (nm) max                     142 non-null    float64
 18  Selectivity (% N_)                         142 non-null    float64
 19  Overall Flow Rate (mL/min)                 142 non-null    object 
 20  Catalyst Weight (g)                        142 non-null    float64
 21  NO Concentration (ppm)                     142 non-null    object 
 22  HC Concentration (ppm)                     142 non-null    object 
 23  Pore Structure                             142 non-null    object 
 24  Acidity                                    142 non-null    object 
 25  Key Reactants                              142 non-null    object 
 26  Temperature Stability (K)                  142 non-null    float64
 27  Deactivation Resistance                    142 non-null    object 
 28  Electron Transfer Properties               142 non-null    object 
 29  Morphology                                 142 non-null    object 
 30  Temperature Range (K)                      142 non-null    object 
 31  NOx Conversion (%)                         142 non-null    float64
dtypes: float64(17), object(15)
memory usage: 36.6+ KB

#converting object numerical columns to float datatype
df['Temperature Range (K)'] = df['Temperature Range (K)'].str.replace(',', '').astype(float)
df['HC Concentration (ppm)'] = df['HC Concentration (ppm)'].str.replace(',', '').astype(float)
df['NO Concentration (ppm)'] = df['NO Concentration (ppm)'].str.replace(',', '').astype(float)
df['Catalyst Molecular Weight (g/mol)'] = df['Catalyst Molecular Weight (g/mol)'].str.replace(',', '').astype(float)
df['Overall Flow Rate (mL/min)'] = df['Overall Flow Rate (mL/min)'].str.replace(',', '').astype(float)

# Visualizing null values using a heatmap
plt.figure(figsize=(5, 3))  # Set the size of the plot
sns.heatmap(df.isnull(), cmap='viridis', cbar=False, yticklabels=False, xticklabels=True)

# Display the plot
plt.title("Heatmap of Null Values")
plt.show()

# Creating subplots for different categorical relationships
fig, axes = plt.subplots(2, 3, figsize=(20, 10))

# Plotting distributions of Pore Structure against various categorical features
sns.histplot(data=df, x="Pore Structure", hue="Support Type", multiple="stack", ax=axes[0, 0])
sns.histplot(data=df, x="Pore Structure", hue="Support Interaction Level", multiple="stack", ax=axes[0, 1])
sns.histplot(data=df, x="Pore Structure", hue="Interaction Type", multiple="stack", ax=axes[1, 0])
sns.histplot(data=df, x="Pore Structure", hue="Acidity", multiple="stack", ax=axes[1, 1])
sns.histplot(data=df, x="Pore Structure", hue="Key Reactants", multiple="stack", ax=axes[0, 2])
sns.histplot(data=df, x="Pore Structure", hue="Morphology", multiple="stack", ax=axes[1, 2])

plt.tight_layout()  # Adjust layout to prevent overlap
plt.show()

# Repeat the above plots for 'Morphology' feature
fig, axes = plt.subplots(2, 3, figsize=(20, 10))

sns.histplot(data=df, x="Morphology", hue="Support Type", multiple="stack", ax=axes[0, 0])
sns.histplot(data=df, x="Morphology", hue="Support Interaction Level", multiple="stack", ax=axes[0, 1])
sns.histplot(data=df, x="Morphology", hue="Interaction Type", multiple="stack", ax=axes[1, 0])
sns.histplot(data=df, x="Morphology", hue="Acidity", multiple="stack", ax=axes[1, 1])
sns.histplot(data=df, x="Morphology", hue="Key Reactants", multiple="stack", ax=axes[0, 2])
sns.histplot(data=df, x="Morphology", hue="Electron Transfer Properties", multiple="stack", ax=axes[1, 2])

plt.tight_layout()  # Adjust layout to prevent overlap
plt.show()

from sklearn.preprocessing import LabelEncoder

# Create a copy of the dataframe for encoding
df2 = df.copy()

# Initialize LabelEncoder
label_encoder = LabelEncoder()

# List of categorical columns to encode
categorical_columns = [
    'Interaction Type',
    'Support Interaction Level',
    'Support Type',
    'Pore Structure',
    'Acidity',
    'Key Reactants',
    'Deactivation Resistance',
    'Electron Transfer Properties',
    'Morphology'
]

# Apply LabelEncoder to each categorical column
for column in categorical_columns:
    df2[column] = label_encoder.fit_transform(df2[column])

# Display the first few rows of the encoded dataset
df2.head()

from mpl_toolkits.mplot3d import Axes3D

# 3D Plot: Surface Area, Pore Volume, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Surface Area (m_/g) min'], 
    df2['Pore Volume (cm_/g) min'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='viridis'
)

ax.set_xlabel('Surface Area (m²/g) min')
ax.set_ylabel('Pore Volume (cm³/g) min')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Surface Area, Pore Volume, and NOx Conversion")
plt.show()

# 3D Plot: Optimal Temperature, Temperature Stability, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Optimal Temperature (K) min'], 
    df2['Temperature Stability (K)'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='plasma'
)

ax.set_xlabel('Optimal Temperature (K) min')
ax.set_ylabel('Temperature Stability (K)')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Optimal Temperature, Temperature Stability, and NOx Conversion")
plt.show()

# 3D Plot: Catalyst Electronegativity, Support Electronegativity, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Catalyst Electronegativity'], 
    df2['Support Material Electronegativity'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='cool'
)

ax.set_xlabel('Catalyst Electronegativity')
ax.set_ylabel('Support Material Electronegativity')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Catalyst and Support Electronegativity and NOx Conversion")
plt.show()

# 3D Plot: Catalyst Molecular Weight, Electronegativity, and Temperature Range
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Catalyst Molecular Weight (g/mol)'], 
    df2['Catalyst Electronegativity'], 
    df2['Temperature Range (K)'], 
    c=df2['Temperature Range (K)'], cmap='coolwarm'
)

ax.set_xlabel('Catalyst Molecular Weight (g/mol)')
ax.set_ylabel('Catalyst Electronegativity')
ax.set_zlabel('Temperature Range (K)')
ax.set_title("3D Plot: Catalyst Molecular Weight, Electronegativity, and Temperature Range")
plt.show()

# 3D Plot: Surface Area, Particle Size, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Surface Area (m_/g) min'], 
    df2['Particle Size (nm) min'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='viridis'
)

ax.set_xlabel('Surface Area (m²/g) min')
ax.set_ylabel('Particle Size (nm) min')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Surface Area, Particle Size, and NOx Conversion")
plt.show()

# 3D Plot: Acidity, Reductive Level, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Acidity'], 
    df2['Reductive Level'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='plasma'
)

ax.set_xlabel('Acidity')
ax.set_ylabel('Reductive Level')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Acidity, Reductive Level, and NOx Conversion")
plt.show()

# 3D Plot: Key Reactants, Flow Rate, and NOx Conversion
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(
    df2['Key Reactants'], 
    df2['Overall Flow Rate (mL/min)'], 
    df2['NOx Conversion (%)'], 
    c=df2['NOx Conversion (%)'], cmap='inferno'
)

ax.set_xlabel('Key Reactants')
ax.set_ylabel('Overall Flow Rate (mL/min)')
ax.set_zlabel('NOx Conversion (%)')
ax.set_title("3D Plot: Key Reactants, Flow Rate, and NOx Conversion")
plt.show()

# Dropping the 'Catalyst Used' column from the dataset
df2 = df2.drop(['Catalyst Used'], axis=1)

# Displaying the first few rows of the updated dataset to confirm the drop
df2.head()

# Correlation matrix
corr_matrix = df2.corr()

# Create a heatmap to visualize the correlation matrix
plt.figure(figsize=(15, 15))
sns.heatmap(corr_matrix, cmap='coolwarm', fmt='.2f', cbar=True, square=True)

# Show the plot
plt.title('Correlation Matrix')
plt.show()

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# Initialize the StandardScaler
scaler = StandardScaler()

# Standardize the data
df_scaled = scaler.fit_transform(df2)

# Apply PCA
pca = PCA()
pca.fit(df_scaled)

# Explained variances (individual and cumulative)
explained_variance = pca.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance)

# Plot individual and cumulative explained variance
plt.figure(figsize=(10, 6))

# Plot individual variance
plt.bar(range(1, len(explained_variance) + 1), explained_variance, alpha=0.6, label='Individual Explained Variance')

# Create a secondary y-axis to plot the cumulative variance
ax2 = plt.gca().twinx()
ax2.plot(range(1, len(cumulative_variance) + 1), cumulative_variance, marker='o', color='b', label='Cumulative Explained Variance')

# Labels and title
plt.xlabel('Principal Component')
plt.ylabel('Individual Variance Explained')
ax2.set_ylabel('Cumulative Variance Explained')
plt.title('PCA: Individual and Cumulative Variance Explained')

# Add legends
plt.legend(loc='upper left')
ax2.legend(loc='upper right')

# Show the plot
plt.tight_layout()
plt.show()

# Apply PCA with first 10 components
pca = PCA(n_components=10)
pca.fit(df_scaled)

# Get the first 10 principal components
pca_components = pca.transform(df_scaled)[:, :10]

# Create a DataFrame for the principal components
pc_df = pd.DataFrame(pca_components, columns=[f'PC{i+1}' for i in range(10)])

# Create the correlation matrix between the original features and the first 10 principal components
correlation_matrix = pd.DataFrame(
    np.corrcoef(df_scaled.T, pc_df.T, rowvar=True)[:df_scaled.shape[1], df_scaled.shape[1]:],
    index=df2.columns,
    columns=pc_df.columns
)

# Plot the correlation matrix using a heatmap
plt.figure(figsize=(10, 10))
sns.heatmap(correlation_matrix, cmap='coolwarm', fmt='.2f', cbar=True, square=True)

# Add title and labels
plt.title('Correlation Matrix: Original Features vs First 10 Principal Components')
plt.show()

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

# Split the data into features and targets
X = df2.iloc[:, :29]
y = df2.iloc[:, 29:31]

# Split into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# Display the shapes of the splits
print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)

(113, 29) (29, 29) (113, 2) (29, 2)

# Initialize the Linear Regression model
model = LinearRegression()

# Train for Temperature Range (K)
model.fit(X_train, y_train['Temperature Range (K)'])
y_pred_1_train = model.predict(X_train)  # Prediction on training data
y_pred_1_test = model.predict(X_test)    # Prediction on testing data

# Train for NOx Conversion (%)
model.fit(X_train, y_train['NOx Conversion (%)'])
y_pred_2_train = model.predict(X_train)  # Prediction on training data
y_pred_2_test = model.predict(X_test)    # Prediction on testing data

# Calculate MSE and R-squared for Temperature Range
mse_1_train = mean_squared_error(y_train['Temperature Range (K)'], y_pred_1_train)
r2_1_train = r2_score(y_train['Temperature Range (K)'], y_pred_1_train)
mse_1_test = mean_squared_error(y_test['Temperature Range (K)'], y_pred_1_test)
r2_1_test = r2_score(y_test['Temperature Range (K)'], y_pred_1_test)

# Calculate MSE and R-squared for NOx Conversion
mse_2_train = mean_squared_error(y_train['NOx Conversion (%)'], y_pred_2_train)
r2_2_train = r2_score(y_train['NOx Conversion (%)'], y_pred_2_train)
mse_2_test = mean_squared_error(y_test['NOx Conversion (%)'], y_pred_2_test)
r2_2_test = r2_score(y_test['NOx Conversion (%)'], y_pred_2_test)

# Add performance metrics to the dictionary
addtodf('LinearRegression_train_NOx', mse_2_train, r2_2_train)
addtodf('LinearRegression_train_temp', mse_1_train, r2_1_train)
addtodf('LinearRegression_test_NOx', mse_2_test, r2_2_test)
addtodf('LinearRegression_test_temp', mse_1_test, r2_1_test)

# Plotting predictions
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Temperature Range (K)
sns.scatterplot(x=y_train['Temperature Range (K)'], y=y_pred_1_train, ax=axes[0], color='blue', label='Train Data')
sns.scatterplot(x=y_test['Temperature Range (K)'], y=y_pred_1_test, ax=axes[0], color='red', label='Test Data')
axes[0].plot(
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    'k--', lw=2
)
axes[0].set_title('Temperature Range (K) - Predicted vs Actual')
axes[0].set_xlabel('Actual Temperature Range (K)')
axes[0].set_ylabel('Predicted Temperature Range (K)')
axes[0].legend()

# NOx Conversion (%)
sns.scatterplot(x=y_train['NOx Conversion (%)'], y=y_pred_2_train, ax=axes[1], color='blue', label='Train Data')
sns.scatterplot(x=y_test['NOx Conversion (%)'], y=y_pred_2_test, ax=axes[1], color='red', label='Test Data')
axes[1].plot(
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    'k--', lw=2
)
axes[1].set_title('NOx Conversion (%) - Predicted vs Actual')
axes[1].set_xlabel('Actual NOx Conversion (%)')
axes[1].set_ylabel('Predicted NOx Conversion (%)')
axes[1].legend()

plt.tight_layout()
plt.show()

# Output performance metrics
print("For Temperature Range (K):")
print(f"Training MSE: {mse_1_train}")
print(f"Training R-squared: {r2_1_train}")
print(f"Test MSE: {mse_1_test}")
print(f"Test R-squared: {r2_1_test}")

print("\nFor NOx Conversion (%):")
print(f"Training MSE: {mse_2_train}")
print(f"Training R-squared: {r2_2_train}")
print(f"Test MSE: {mse_2_test}")
print(f"Test R-squared: {r2_2_test}")

For Temperature Range (K):
Training MSE: 16177.587784781055
Training R-squared: 0.5407359206687252
Test MSE: 45399.76260510736
Test R-squared: -0.021614175796591528

For NOx Conversion (%):
Training MSE: 489.3767142516031
Training R-squared: 0.23540240102346188
Test MSE: 1052.375345293837
Test R-squared: -1.2801781045029905

from sklearn.ensemble import RandomForestRegressor

# Initialize the RandomForestRegressor
model = RandomForestRegressor(random_state=12)

# Train for Temperature Range (K)
model.fit(X_train, y_train['Temperature Range (K)'])
y_pred_1_train = model.predict(X_train)  # Prediction on training data
y_pred_1_test = model.predict(X_test)    # Prediction on testing data

# Store feature importance for Temperature Range (K)
importance_1 = model.feature_importances_

# Train for NOx Conversion (%)
model.fit(X_train, y_train['NOx Conversion (%)'])
y_pred_2_train = model.predict(X_train)  # Prediction on training data
y_pred_2_test = model.predict(X_test)    # Prediction on testing data

# Store feature importance for NOx Conversion (%)
importance_2 = model.feature_importances_

# Calculate MSE and R-squared for Temperature Range
mse_1_train = mean_squared_error(y_train['Temperature Range (K)'], y_pred_1_train)
r2_1_train = r2_score(y_train['Temperature Range (K)'], y_pred_1_train)
mse_1_test = mean_squared_error(y_test['Temperature Range (K)'], y_pred_1_test)
r2_1_test = r2_score(y_test['Temperature Range (K)'], y_pred_1_test)

# Calculate MSE and R-squared for NOx Conversion
mse_2_train = mean_squared_error(y_train['NOx Conversion (%)'], y_pred_2_train)
r2_2_train = r2_score(y_train['NOx Conversion (%)'], y_pred_2_train)
mse_2_test = mean_squared_error(y_test['NOx Conversion (%)'], y_pred_2_test)
r2_2_test = r2_score(y_test['NOx Conversion (%)'], y_pred_2_test)

# Add performance metrics to the dictionary
addtodf('RandomForestRegressor_train_NOx', mse_2_train, r2_2_train)
addtodf('RandomForestRegressor_train_temp', mse_1_train, r2_1_train)
addtodf('RandomForestRegressor_test_NOx', mse_2_test, r2_2_test)
addtodf('RandomForestRegressor_test_temp', mse_1_test, r2_1_test)

# Plotting predictions
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Temperature Range (K)
sns.scatterplot(x=y_train['Temperature Range (K)'], y=y_pred_1_train, ax=axes[0], color='blue', label='Train Data')
sns.scatterplot(x=y_test['Temperature Range (K)'], y=y_pred_1_test, ax=axes[0], color='red', label='Test Data')
axes[0].plot(
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    'k--', lw=2
)
axes[0].set_title('Temperature Range (K) - Predicted vs Actual')
axes[0].set_xlabel('Actual Temperature Range (K)')
axes[0].set_ylabel('Predicted Temperature Range (K)')
axes[0].legend()

# NOx Conversion (%)
sns.scatterplot(x=y_train['NOx Conversion (%)'], y=y_pred_2_train, ax=axes[1], color='blue', label='Train Data')
sns.scatterplot(x=y_test['NOx Conversion (%)'], y=y_pred_2_test, ax=axes[1], color='red', label='Test Data')
axes[1].plot(
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    'k--', lw=2
)
axes[1].set_title('NOx Conversion (%) - Predicted vs Actual')
axes[1].set_xlabel('Actual NOx Conversion (%)')
axes[1].set_ylabel('Predicted NOx Conversion (%)')
axes[1].legend()

plt.tight_layout()
plt.show()

# Output performance metrics
print("For Temperature Range (K):")
print(f"Training MSE: {mse_1_train}")
print(f"Training R-squared: {r2_1_train}")
print(f"Test MSE: {mse_1_test}")
print(f"Test R-squared: {r2_1_test}")

print("\nFor NOx Conversion (%):")
print(f"Training MSE: {mse_2_train}")
print(f"Training R-squared: {r2_2_train}")
print(f"Test MSE: {mse_2_test}")
print(f"Test R-squared: {r2_2_test}")

For Temperature Range (K):
Training MSE: 2784.5424583317754
Training R-squared: 0.9209498755007439
Test MSE: 19385.292694572596
Test R-squared: 0.5637798816041928

For NOx Conversion (%):
Training MSE: 111.20186696779744
Training R-squared: 0.8262592436272467
Test MSE: 463.76222983716474
Test R-squared: -0.004832056261190321

# Plotting Feature Importance for Temperature Range (K)
plt.figure(figsize=(12, 8))
sns.barplot(x=importance_1, y=X_train.columns, palette="viridis")
plt.title("Feature Importance for Temperature Range (K)")
plt.xlabel("Importance")
plt.ylabel("Features")
plt.show()

# Plotting Feature Importance for NOx Conversion (%)
plt.figure(figsize=(12, 8))
sns.barplot(x=importance_2, y=X_train.columns, palette="plasma")
plt.title("Feature Importance for NOx Conversion (%)")
plt.xlabel("Importance")
plt.ylabel("Features")
plt.show()

import xgboost as xgb

# Initialize the XGBRegressor
model = xgb.XGBRegressor(random_state=12, booster='gbtree', max_depth=12, learning_rate=0.01, n_estimators=200)

# Train for Temperature Range (K)
model.fit(X_train, y_train['Temperature Range (K)'])
y_pred_1_train = model.predict(X_train)  # Prediction on training data
y_pred_1_test = model.predict(X_test)    # Prediction on testing data

# Store feature importance for Temperature Range (K)
importance_1 = model.get_booster().get_score(importance_type='weight')  # 'weight', 'gain', or 'cover'

# Train for NOx Conversion (%)
model.fit(X_train, y_train['NOx Conversion (%)'])
y_pred_2_train = model.predict(X_train)  # Prediction on training data
y_pred_2_test = model.predict(X_test)    # Prediction on testing data

# Store feature importance for NOx Conversion (%)
importance_2 = model.get_booster().get_score(importance_type='weight')  # 'weight', 'gain', or 'cover'

# Calculate MSE and R-squared for Temperature Range
mse_1_train = mean_squared_error(y_train['Temperature Range (K)'], y_pred_1_train)
r2_1_train = r2_score(y_train['Temperature Range (K)'], y_pred_1_train)
mse_1_test = mean_squared_error(y_test['Temperature Range (K)'], y_pred_1_test)
r2_1_test = r2_score(y_test['Temperature Range (K)'], y_pred_1_test)

# Calculate MSE and R-squared for NOx Conversion
mse_2_train = mean_squared_error(y_train['NOx Conversion (%)'], y_pred_2_train)
r2_2_train = r2_score(y_train['NOx Conversion (%)'], y_pred_2_train)
mse_2_test = mean_squared_error(y_test['NOx Conversion (%)'], y_pred_2_test)
r2_2_test = r2_score(y_test['NOx Conversion (%)'], y_pred_2_test)

# Add performance metrics to the dictionary
addtodf('XGBRegressor_train_NOx', mse_2_train, r2_2_train)
addtodf('XGBRegressor_train_temp', mse_1_train, r2_1_train)
addtodf('XGBRegressor_test_NOx', mse_2_test, r2_2_test)
addtodf('XGBRegressor_test_temp', mse_1_test, r2_1_test)

# Plotting predictions
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Temperature Range (K)
sns.scatterplot(x=y_train['Temperature Range (K)'], y=y_pred_1_train, ax=axes[0], color='blue', label='Train Data')
sns.scatterplot(x=y_test['Temperature Range (K)'], y=y_pred_1_test, ax=axes[0], color='red', label='Test Data')
axes[0].plot(
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    [y_train['Temperature Range (K)'].min(), y_train['Temperature Range (K)'].max()],
    'k--', lw=2
)
axes[0].set_title('Temperature Range (K) - Predicted vs Actual')
axes[0].set_xlabel('Actual Temperature Range (K)')
axes[0].set_ylabel('Predicted Temperature Range (K)')
axes[0].legend()

# NOx Conversion (%)
sns.scatterplot(x=y_train['NOx Conversion (%)'], y=y_pred_2_train, ax=axes[1], color='blue', label='Train Data')
sns.scatterplot(x=y_test['NOx Conversion (%)'], y=y_pred_2_test, ax=axes[1], color='red', label='Test Data')
axes[1].plot(
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    [y_train['NOx Conversion (%)'].min(), y_train['NOx Conversion (%)'].max()],
    'k--', lw=2
)
axes[1].set_title('NOx Conversion (%) - Predicted vs Actual')
axes[1].set_xlabel('Actual NOx Conversion (%)')
axes[1].set_ylabel('Predicted NOx Conversion (%)')
axes[1].legend()

plt.tight_layout()
plt.show()

# Output performance metrics
print("For Temperature Range (K):")
print(f"Training MSE: {mse_1_train}")
print(f"Training R-squared: {r2_1_train}")
print(f"Test MSE: {mse_1_test}")
print(f"Test R-squared: {r2_1_test}")

print("\nFor NOx Conversion (%):")
print(f"Training MSE: {mse_2_train}")
print(f"Training R-squared: {r2_2_train}")
print(f"Test MSE: {mse_2_test}")
print(f"Test R-squared: {r2_2_test}")

For Temperature Range (K):
Training MSE: 2101.979284587625
Training R-squared: 0.9403270998277193
Test MSE: 25962.91766276785
Test R-squared: 0.41576600388788265

For NOx Conversion (%):
Training MSE: 58.873610278569494
Training R-squared: 0.9080164222138871
Test MSE: 534.7799242252285
Test R-squared: -0.15870585471162268

# Plotting Feature Importance for Temperature Range (K)
plt.figure(figsize=(12, 8))
sns.barplot(x=list(importance_1.values()), y=list(importance_1.keys()), palette="viridis")
plt.title("Feature Importance for Temperature Range (K)")
plt.xlabel("Importance (Weight)")
plt.ylabel("Features")
plt.show()

# Plotting Feature Importance for NOx Conversion (%)
plt.figure(figsize=(12, 8))
sns.barplot(x=list(importance_2.values()), y=list(importance_2.keys()), palette="plasma")
plt.title("Feature Importance for NOx Conversion (%)")
plt.xlabel("Importance (Weight)")
plt.ylabel("Features")
plt.show()

from pyDOE import lhs
from mpl_toolkits.mplot3d import Axes3D

# Define feature matrix and target variables
X = df2.iloc[:, :29]
y_temp = df2['Temperature Range (K)']
y_nox = df2['NOx Conversion (%)']

# Initialize lists to store predictions from each model
all_predictions_temp = []
all_predictions_nox = []

# Train multiple XGBoost models with different splits
for i in range(4):
    # Split the data for Temperature Range (K)
    X_train_split, X_test_split, y_train_temp_split, y_test_temp_split = train_test_split(
        X, y_temp, test_size=0.4, random_state=i
    )
    model_temp = xgb.XGBRegressor(random_state=12)
    model_temp.fit(X_train_split, y_train_temp_split)
    y_pred_temp = model_temp.predict(X)
    all_predictions_temp.append(y_pred_temp)
    
    # Split the data for NOx Conversion (%)
    _, _, y_train_nox_split, y_test_nox_split = train_test_split(
        X, y_nox, test_size=0.4, random_state=i
    )
    model_nox = xgb.XGBRegressor(random_state=12)
    model_nox.fit(X_train_split, y_train_nox_split)
    y_pred_nox = model_nox.predict(X)
    all_predictions_nox.append(y_pred_nox)

# Convert predictions lists to 2D arrays (n_samples, number_of_models)
all_predictions_temp = np.column_stack(all_predictions_temp)
all_predictions_nox = np.column_stack(all_predictions_nox)

# Calculate the mean of the predictions
mean_predictions_temp = np.mean(all_predictions_temp, axis=1)
mean_predictions_nox = np.mean(all_predictions_nox, axis=1)

# Evaluate the averaged predictions
mse_temp = mean_squared_error(y_temp, mean_predictions_temp)
r2_temp = r2_score(y_temp, mean_predictions_temp)
mse_nox = mean_squared_error(y_nox, mean_predictions_nox)
r2_nox = r2_score(y_nox, mean_predictions_nox)

# Add performance metrics to the dictionary
addtodf('XGBRegressor_100_NOx', mse_nox, r2_nox)
addtodf('XGBRegressor_100_temp', mse_temp, r2_temp)

# Display performance metrics
print("=== Model Performance for Temperature Range (K) ===")
print("Mean Squared Error (MSE):", mse_temp)
print("R-squared (R²):", r2_temp)

print("\n=== Model Performance for NOx Conversion (%) ===")
print("Mean Squared Error (MSE):", mse_nox)
print("R-squared (R²):", r2_nox)

=== Model Performance for Temperature Range (K) ===
Mean Squared Error (MSE): 4562.457067962988
R-squared (R²): 0.8774906370383956

=== Model Performance for NOx Conversion (%) ===
Mean Squared Error (MSE): 160.01272822998376
R-squared (R²): 0.7355008392408131

# Plot the averaged predictions vs actual values for both targets
plt.figure(figsize=(14, 6))

# Temperature Range (K)
plt.subplot(1, 2, 1)
plt.scatter(y_temp, mean_predictions_temp, color='blue', label='Averaged Predictions')
plt.plot(
    [y_temp.min(), y_temp.max()],
    [y_temp.min(), y_temp.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual Temperature Range (K)")
plt.ylabel("Predicted Temperature Range (K)")
plt.title("Averaged Model Predictions vs Actual Values for Temperature Range (K)")
plt.legend()
plt.grid(True)

# NOx Conversion (%)
plt.subplot(1, 2, 2)
plt.scatter(y_nox, mean_predictions_nox, color='green', label='Averaged Predictions')
plt.plot(
    [y_nox.min(), y_nox.max()],
    [y_nox.min(), y_nox.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual NOx Conversion (%)")
plt.ylabel("Predicted NOx Conversion (%)")
plt.title("Averaged Model Predictions vs Actual Values for NOx Conversion (%)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Initialize lists to store predictions from each model
all_predictions_temp = []
all_predictions_nox = []

# Train multiple Random Forest models with different splits
for i in range(4):
    # Split the data for Temperature Range (K)
    X_train_split, X_test_split, y_train_temp_split, y_test_temp_split = train_test_split(
        X, y_temp, test_size=0.4, random_state=i
    )
    model_temp = RandomForestRegressor(random_state=12)
    model_temp.fit(X_train_split, y_train_temp_split)
    y_pred_temp = model_temp.predict(X)
    all_predictions_temp.append(y_pred_temp)
    
    # Split the data for NOx Conversion (%)
    _, _, y_train_nox_split, y_test_nox_split = train_test_split(
        X, y_nox, test_size=0.4, random_state=i
    )
    model_nox = RandomForestRegressor(random_state=12)
    model_nox.fit(X_train_split, y_train_nox_split)
    y_pred_nox = model_nox.predict(X)
    all_predictions_nox.append(y_pred_nox)

# Convert predictions lists to 2D arrays (n_samples, number_of_models)
all_predictions_temp = np.column_stack(all_predictions_temp)
all_predictions_nox = np.column_stack(all_predictions_nox)

# Calculate the mean of the predictions
mean_predictions_temp = np.mean(all_predictions_temp, axis=1)
mean_predictions_nox = np.mean(all_predictions_nox, axis=1)

# Evaluate the averaged predictions
mse_temp = mean_squared_error(y_temp, mean_predictions_temp)
r2_temp = r2_score(y_temp, mean_predictions_temp)
mse_nox = mean_squared_error(y_nox, mean_predictions_nox)
r2_nox = r2_score(y_nox, mean_predictions_nox)

# Add performance metrics to the dictionary
addtodf('RandomForestRegressor_100_NOx', mse_nox, r2_nox)
addtodf('RandomForestRegressor_100_temp', mse_temp, r2_temp)

# Display performance metrics
print("=== Model Performance for Temperature Range (K) ===")
print("Mean Squared Error (MSE):", mse_temp)
print("R-squared (R²):", r2_temp)

print("\n=== Model Performance for NOx Conversion (%) ===")
print("Mean Squared Error (MSE):", mse_nox)
print("R-squared (R²):", r2_nox)

=== Model Performance for Temperature Range (K) ===
Mean Squared Error (MSE): 7479.103760679326
R-squared (R²): 0.7991739486868971

=== Model Performance for NOx Conversion (%) ===
Mean Squared Error (MSE): 251.4317256501125
R-squared (R²): 0.5843863100246347

# Plot the averaged predictions vs actual values for both targets
plt.figure(figsize=(14, 6))

# Temperature Range (K)
plt.subplot(1, 2, 1)
plt.scatter(y_temp, mean_predictions_temp, color='blue', label='Averaged Predictions')
plt.plot(
    [y_temp.min(), y_temp.max()],
    [y_temp.min(), y_temp.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual Temperature Range (K)")
plt.ylabel("Predicted Temperature Range (K)")
plt.title("Averaged Model Predictions vs Actual Values for Temperature Range (K)")
plt.legend()
plt.grid(True)

# NOx Conversion (%)
plt.subplot(1, 2, 2)
plt.scatter(y_nox, mean_predictions_nox, color='green', label='Averaged Predictions')
plt.plot(
    [y_nox.min(), y_nox.max()],
    [y_nox.min(), y_nox.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual NOx Conversion (%)")
plt.ylabel("Predicted NOx Conversion (%)")
plt.title("Averaged Model Predictions vs Actual Values for NOx Conversion (%)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Initialize lists to store predictions from each model
all_predictions_temp = []
all_predictions_nox = []

# Train multiple Linear Regression models with different splits
for i in range(4):
    # Split the data for Temperature Range (K)
    X_train_split, X_test_split, y_train_temp_split, y_test_temp_split = train_test_split(
        X, y_temp, test_size=0.4, random_state=i
    )
    model_temp = LinearRegression()
    model_temp.fit(X_train_split, y_train_temp_split)
    y_pred_temp = model_temp.predict(X)
    all_predictions_temp.append(y_pred_temp)
    
    # Split the data for NOx Conversion (%)
    _, _, y_train_nox_split, y_test_nox_split = train_test_split(
        X, y_nox, test_size=0.4, random_state=i
    )
    model_nox = LinearRegression()
    model_nox.fit(X_train_split, y_train_nox_split)
    y_pred_nox = model_nox.predict(X)
    all_predictions_nox.append(y_pred_nox)

# Convert predictions lists to 2D arrays (n_samples, number_of_models)
all_predictions_temp = np.column_stack(all_predictions_temp)
all_predictions_nox = np.column_stack(all_predictions_nox)

# Calculate the mean of the predictions
mean_predictions_temp = np.mean(all_predictions_temp, axis=1)
mean_predictions_nox = np.mean(all_predictions_nox, axis=1)

# Evaluate the averaged predictions
mse_temp = mean_squared_error(y_temp, mean_predictions_temp)
r2_temp = r2_score(y_temp, mean_predictions_temp)
mse_nox = mean_squared_error(y_nox, mean_predictions_nox)
r2_nox = r2_score(y_nox, mean_predictions_nox)

# Add performance metrics to the dictionary
addtodf('LinearRegression_100_NOx', mse_nox, r2_nox)
addtodf('LinearRegression_100_temp', mse_temp, r2_temp)

# Display performance metrics
print("=== Model Performance for Temperature Range (K) ===")
print("Mean Squared Error (MSE):", mse_temp)
print("R-squared (R²):", r2_temp)

print("\n=== Model Performance for NOx Conversion (%) ===")
print("Mean Squared Error (MSE):", mse_nox)
print("R-squared (R²):", r2_nox)

=== Model Performance for Temperature Range (K) ===
Mean Squared Error (MSE): 20517.876858319272
R-squared (R²): 0.44906176961362454

=== Model Performance for NOx Conversion (%) ===
Mean Squared Error (MSE): 553.0334189410929
R-squared (R²): 0.08584225267716572

# Plot the averaged predictions vs actual values for both targets
plt.figure(figsize=(14, 6))

# Temperature Range (K)
plt.subplot(1, 2, 1)
plt.scatter(y_temp, mean_predictions_temp, color='blue', label='Averaged Predictions')
plt.plot(
    [y_temp.min(), y_temp.max()],
    [y_temp.min(), y_temp.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual Temperature Range (K)")
plt.ylabel("Predicted Temperature Range (K)")
plt.title("Averaged Model Predictions vs Actual Values for Temperature Range (K)")
plt.legend()
plt.grid(True)

# NOx Conversion (%)
plt.subplot(1, 2, 2)
plt.scatter(y_nox, mean_predictions_nox, color='green', label='Averaged Predictions')
plt.plot(
    [y_nox.min(), y_nox.max()],
    [y_nox.min(), y_nox.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual NOx Conversion (%)")
plt.ylabel("Predicted NOx Conversion (%)")
plt.title("Averaged Model Predictions vs Actual Values for NOx Conversion (%)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

from sklearn.ensemble import GradientBoostingRegressor

# Initialize lists to store predictions from each model
all_predictions_temp = []
all_predictions_nox = []

# Train multiple Gradient Boosting models with different splits
for i in range(4):
    # Split the data for Temperature Range (K)
    X_train_split, X_test_split, y_train_temp_split, y_test_temp_split = train_test_split(
        X, y_temp, test_size=0.4, random_state=i
    )
    model_temp = GradientBoostingRegressor(random_state=12)
    model_temp.fit(X_train_split, y_train_temp_split)
    y_pred_temp = model_temp.predict(X)
    all_predictions_temp.append(y_pred_temp)
    
    # Split the data for NOx Conversion (%)
    _, _, y_train_nox_split, y_test_nox_split = train_test_split(
        X, y_nox, test_size=0.4, random_state=i
    )
    model_nox = GradientBoostingRegressor(random_state=12)
    model_nox.fit(X_train_split, y_train_nox_split)
    y_pred_nox = model_nox.predict(X)
    all_predictions_nox.append(y_pred_nox)

# Convert predictions lists to 2D arrays (n_samples, number_of_models)
all_predictions_temp = np.column_stack(all_predictions_temp)
all_predictions_nox = np.column_stack(all_predictions_nox)

# Calculate the mean of the predictions
mean_predictions_temp = np.mean(all_predictions_temp, axis=1)
mean_predictions_nox = np.mean(all_predictions_nox, axis=1)

# Evaluate the averaged predictions
mse_temp = mean_squared_error(y_temp, mean_predictions_temp)
r2_temp = r2_score(y_temp, mean_predictions_temp)
mse_nox = mean_squared_error(y_nox, mean_predictions_nox)
r2_nox = r2_score(y_nox, mean_predictions_nox)

# Add performance metrics to the dictionary
addtodf('GradientBoostingRegressor_100_NOx', mse_nox, r2_nox)
addtodf('GradientBoostingRegressor_100_temp', mse_temp, r2_temp)

# Display performance metrics
print("=== Model Performance for Temperature Range (K) ===")
print("Mean Squared Error (MSE):", mse_temp)
print("R-squared (R²):", r2_temp)

print("\n=== Model Performance for NOx Conversion (%) ===")
print("Mean Squared Error (MSE):", mse_nox)
print("R-squared (R²):", r2_nox)

=== Model Performance for Temperature Range (K) ===
Mean Squared Error (MSE): 4825.3198456269265
R-squared (R²): 0.8704323456488604

=== Model Performance for NOx Conversion (%) ===
Mean Squared Error (MSE): 179.0140044239983
R-squared (R²): 0.7040919528211853

# Plot the averaged predictions vs actual values for both targets
plt.figure(figsize=(14, 6))

# Temperature Range (K)
plt.subplot(1, 2, 1)
plt.scatter(y_temp, mean_predictions_temp, color='blue', label='Averaged Predictions')
plt.plot(
    [y_temp.min(), y_temp.max()],
    [y_temp.min(), y_temp.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual Temperature Range (K)")
plt.ylabel("Predicted Temperature Range (K)")
plt.title("Averaged Model Predictions vs Actual Values for Temperature Range (K)")
plt.legend()
plt.grid(True)

# NOx Conversion (%)
plt.subplot(1, 2, 2)
plt.scatter(y_nox, mean_predictions_nox, color='green', label='Averaged Predictions')
plt.plot(
    [y_nox.min(), y_nox.max()],
    [y_nox.min(), y_nox.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual NOx Conversion (%)")
plt.ylabel("Predicted NOx Conversion (%)")
plt.title("Averaged Model Predictions vs Actual Values for NOx Conversion (%)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import catboost as cb

# Initialize lists to store predictions from each model
all_predictions_temp = []
all_predictions_nox = []

# Train multiple CatBoost models with different splits
for i in range(100):
    # Split the data for Temperature Range (K)
    X_train_split, X_test_split, y_train_temp_split, y_test_temp_split = train_test_split(
        X, y_temp, test_size=0.4, random_state=i
    )
    model_temp = cb.CatBoostRegressor(silent=True, random_state=12)
    model_temp.fit(X_train_split, y_train_temp_split)
    y_pred_temp = model_temp.predict(X)
    all_predictions_temp.append(y_pred_temp)
    
    # Split the data for NOx Conversion (%)
    _, _, y_train_nox_split, y_test_nox_split = train_test_split(
        X, y_nox, test_size=0.4, random_state=i
    )
    model_nox = cb.CatBoostRegressor(silent=True, random_state=12)
    model_nox.fit(X_train_split, y_train_nox_split)
    y_pred_nox = model_nox.predict(X)
    all_predictions_nox.append(y_pred_nox)

# Convert predictions lists to 2D arrays (n_samples, number_of_models)
all_predictions_temp = np.column_stack(all_predictions_temp)
all_predictions_nox = np.column_stack(all_predictions_nox)

# Calculate the mean of the predictions
mean_predictions_temp = np.mean(all_predictions_temp, axis=1)
mean_predictions_nox = np.mean(all_predictions_nox, axis=1)

# Evaluate the averaged predictions
mse_temp = mean_squared_error(y_temp, mean_predictions_temp)
r2_temp = r2_score(y_temp, mean_predictions_temp)
mse_nox = mean_squared_error(y_nox, mean_predictions_nox)
r2_nox = r2_score(y_nox, mean_predictions_nox)

# Add performance metrics to the dictionary
addtodf('CatBoostRegressor_100_NOx', mse_nox, r2_nox)
addtodf('CatBoostRegressor_100_temp', mse_temp, r2_temp)

# Display performance metrics
print("=== Model Performance for Temperature Range (K) ===")
print("Mean Squared Error (MSE):", mse_temp)
print("R-squared (R²):", r2_temp)

print("\n=== Model Performance for NOx Conversion (%) ===")
print("Mean Squared Error (MSE):", mse_nox)
print("R-squared (R²):", r2_nox)

=== Model Performance for Temperature Range (K) ===
Mean Squared Error (MSE): 3061.498261868303
R-squared (R²): 0.9177938123729019

=== Model Performance for NOx Conversion (%) ===
Mean Squared Error (MSE): 107.94338454475486
R-squared (R²): 0.8215708529101633

# Plot the averaged predictions vs actual values for both targets
plt.figure(figsize=(14, 6))

# Temperature Range (K)
plt.subplot(1, 2, 1)
plt.scatter(y_temp, mean_predictions_temp, color='blue', label='Averaged Predictions')
plt.plot(
    [y_temp.min(), y_temp.max()],
    [y_temp.min(), y_temp.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual Temperature Range (K)")
plt.ylabel("Predicted Temperature Range (K)")
plt.title("Averaged Model Predictions vs Actual Values for Temperature Range (K)")
plt.legend()
plt.grid(True)

# NOx Conversion (%)
plt.subplot(1, 2, 2)
plt.scatter(y_nox, mean_predictions_nox, color='green', label='Averaged Predictions')
plt.plot(
    [y_nox.min(), y_nox.max()],
    [y_nox.min(), y_nox.max()],
    color='red', linewidth=2, label='Ideal Fit'
)
plt.xlabel("Actual NOx Conversion (%)")
plt.ylabel("Predicted NOx Conversion (%)")
plt.title("Averaged Model Predictions vs Actual Values for NOx Conversion (%)")
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Create a DataFrame to display model performance metrics
performance_df = pd.DataFrame(modeldict)
performance_df

from pyDOE import lhs
from mpl_toolkits.mplot3d import Axes3D

# Select a specific catalyst by its ID
catalyst_id = 44

# Extract fixed variable values from the dataset
a1 = df2.iloc[catalyst_id, :].to_dict()['NO Concentration (ppm)']
a2 = df2.iloc[catalyst_id, :].to_dict()['HC Concentration (ppm)']
a3 = df2.iloc[catalyst_id, :].to_dict()['Overall Flow Rate (mL/min)']

# Define the ranges for the parameters based on the selected catalyst
num_samples = 2000  # Number of samples to generate
no_range = (a1 - 1000, a1 + 1000)         # Range for NO Concentration (ppm)
hc_range = (a2 - 1000, a2 + 1000)         # Range for HC Concentration (ppm)
flow_rate_range = (a3 - 1000, a3 + 1000) # Range for Overall Flow Rate (mL/min)

# Perform Latin Hypercube Sampling
lhs_samples = lhs(3, samples=num_samples)  # 3 parameters: NO, HC, Flow Rate

# Scale LHS samples to the specified ranges
df1 = pd.DataFrame({
    'NO Concentration (ppm)': (lhs_samples[:, 0] * (no_range[1] - no_range[0]) + no_range[0]).astype(int),
    'HC Concentration (ppm)': (lhs_samples[:, 1] * (hc_range[1] - hc_range[0]) + hc_range[0]).astype(int),
    'Overall Flow Rate (mL/min)': (lhs_samples[:, 2] * (flow_rate_range[1] - flow_rate_range[0]) + flow_rate_range[0]).astype(int)
})

# Add fixed variables to the new DataFrame
fixed_values = df2.iloc[catalyst_id, :].to_dict()

for column, value in fixed_values.items():
    if column not in ['NO Concentration (ppm)', 'HC Concentration (ppm)', 'Overall Flow Rate (mL/min)']:
        df1[column] = value

# Reorder the columns to match the original dataset
columns_order = X_train.columns
df1 = df1[columns_order]

# Predict NOx conversion and temperature using pre-trained models
predicted_NOx = model_nox.predict(df1)
predicted_temp = model_temp.predict(df1)


# Create a new DataFrame for plotting
df_plot = df1[['NO Concentration (ppm)', 'HC Concentration (ppm)']].copy()
df_plot['Predicted NOx conversion (model_NOx)'] = predicted_NOx
df_plot['Predicted temperature (model_temp)'] = predicted_temp

# Get the original value for catalyst_id from df2
original_values = df2.iloc[catalyst_id, 20:22].to_dict()  # Adjust indices as per your dataset
original_NO = original_values['NO Concentration (ppm)']
original_HC = original_values['HC Concentration (ppm)']

# Predict the NOx conversion and temperature for the original values
original_predicted_NOx = df2.iloc[catalyst_id, :].to_dict()['NOx Conversion (%)']
original_predicted_temp = df2.iloc[catalyst_id, :].to_dict()['Temperature Range (K)']

# Create a 3D plot for NOx conversion vs NO and HC Concentrations
fig = plt.figure(figsize=(12, 10))

# 3D Plot for NOx Conversion
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(
    df_plot['NO Concentration (ppm)'], 
    df_plot['HC Concentration (ppm)'], 
    df_plot['Predicted NOx conversion (model_NOx)'], 
    color='orange', alpha=0.6
)
ax1.scatter(
    original_NO, original_HC, original_predicted_NOx, 
    color='red', s=100, label='Original Value', marker='X'
)
ax1.set_xlabel('NO Concentration (ppm)')
ax1.set_ylabel('HC Concentration (ppm)')
ax1.set_zlabel('NOx Conversion (%)')
ax1.set_title('3D Plot: NOx Conversion vs NO and HC Concentration')
ax1.legend()

# 3D Plot for Temperature
ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(
    df_plot['NO Concentration (ppm)'], 
    df_plot['HC Concentration (ppm)'], 
    df_plot['Predicted temperature (model_temp)'], 
    color='green', alpha=0.6
)
ax2.scatter(
    original_NO, original_HC, original_predicted_temp, 
    color='blue', s=100, label='Original Value', marker='X'
)
ax2.set_xlabel('NO Concentration (ppm)')
ax2.set_ylabel('HC Concentration (ppm)')
ax2.set_zlabel('Temperature (°C)')
ax2.set_title('3D Plot: Temperature vs NO and HC Concentration')
ax2.legend()

plt.tight_layout()
plt.show()

# Define the ranges for parameters with finer granularity
no_range_fine = np.linspace(a1 - 1000, a1 + 1000, 100)  # 100 points for NO Concentration
hc_range_fine = np.linspace(a2 - 1000, a2 + 1000, 100)  # 100 points for HC Concentration

# Create a meshgrid for NO and HC concentrations
NO, HC = np.meshgrid(no_range_fine, hc_range_fine)

# Flatten the grid for passing into models
df1_grid = pd.DataFrame({
    'NO Concentration (ppm)': NO.flatten(),
    'HC Concentration (ppm)': HC.flatten()
})

# Add fixed variables to the new DataFrame
for column, value in fixed_values.items():
    if column not in ['NO Concentration (ppm)', 'HC Concentration (ppm)']:
        df1_grid[column] = value

# Reorder the columns to match the original dataset
df1_grid = df1_grid[columns_order]

# Predict NOx conversion and temperature
predicted_NOx_grid = model_nox.predict(df1_grid)
predicted_temp_grid = model_temp.predict(df1_grid)

# Reshape the predicted values back into the shape of the meshgrid
predicted_NOx_grid_reshaped = predicted_NOx_grid.reshape(NO.shape)
predicted_temp_grid_reshaped = predicted_temp_grid.reshape(NO.shape)

# Create a 3D surface plot for NOx conversion
fig = plt.figure(figsize=(15, 10))

# Surface Plot for NOx Conversion
ax1 = fig.add_subplot(121, projection='3d')
surf1 = ax1.plot_surface(NO, HC, predicted_NOx_grid_reshaped, cmap='YlOrRd', edgecolor='none')
fig.colorbar(surf1, ax=ax1)
ax1.set_xlabel('NO Concentration (ppm)')
ax1.set_ylabel('HC Concentration (ppm)')
ax1.set_zlabel('NOx Conversion (%)')
ax1.set_title('Surface Plot: NOx Conversion vs NO and HC Concentration')

# Surface Plot for Temperature
ax2 = fig.add_subplot(122, projection='3d')
surf2 = ax2.plot_surface(NO, HC, predicted_temp_grid_reshaped, cmap='coolwarm', edgecolor='none')
fig.colorbar(surf2, ax=ax2)
ax2.set_xlabel('NO Concentration (ppm)')
ax2.set_ylabel('HC Concentration (ppm)')
ax2.set_zlabel('Temperature (°C)')
ax2.set_title('Surface Plot: Temperature vs NO and HC Concentration')

plt.tight_layout()
plt.show()

# Function to predict NOx Conversion and Temperature Range based on input concentrations and flow rate
def predict_optimal_conditions(no_concentration_ppm, hc_concentration_ppm, flow_rate_ml_min):
    """
    Predicts NOx Conversion (%) and Temperature Range (K) based on input concentrations and flow rate.
    
    Parameters:
    - no_concentration_ppm (int): Concentration of nitric oxide (ppm).
    - hc_concentration_ppm (int): Concentration of hydrocarbons (ppm).
    - flow_rate_ml_min (int): Overall flow rate (mL/min).
    
    Returns:
    - dict: Predicted NOx Conversion (%) and Temperature Range (K).
    """
    # Create a DataFrame with the input values
    input_data = pd.DataFrame({
        'NO Concentration (ppm)': [no_concentration_ppm],
        'HC Concentration (ppm)': [hc_concentration_ppm],
        'Overall Flow Rate (mL/min)': [flow_rate_ml_min]
    })
    
    # Add fixed features from a selected catalyst (assuming catalyst_id = 44)
    catalyst_id = 44
    fixed_values = df2.iloc[catalyst_id, :].to_dict()
    
    for column, value in fixed_values.items():
        if column not in ['NO Concentration (ppm)', 'HC Concentration (ppm)', 'Overall Flow Rate (mL/min)']:
            input_data[column] = value
    
    # Ensure the input_data has the same order of columns as the training data
    input_data = input_data[X_train.columns]
    
    # Predict NOx Conversion using the best-performing model (e.g., XGBoost Stacking)
    predicted_NOx = model_nox.predict(input_data)[0]
    
    # Predict Temperature Range using the best-performing model (e.g., XGBoost Stacking)
    predicted_temp = model_temp.predict(input_data)[0]
    
    return {
        'Predicted NOx Conversion (%)': round(predicted_NOx, 2),
        'Predicted Temperature Range (K)': round(predicted_temp, 2)
    }

# Example usage of the function
# Input values
no_input = 1500        # NO Concentration in ppm
hc_input = 800         # HC Concentration in ppm
flow_rate_input = 500  # Flow Rate in mL/min

# Get predictions
predictions = predict_optimal_conditions(no_input, hc_input, flow_rate_input)

# Display the predictions
print("Optimal Reaction Conditions:")
print(f"NOx Conversion (%): {predictions['Predicted NOx Conversion (%)']}")
print(f"Temperature Range (K): {predictions['Predicted Temperature Range (K)']}")

Optimal Reaction Conditions:
NOx Conversion (%): 69.22
Temperature Range (K): 671.75